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Application of a dynamic subgrid-scale 
model to turbulent recirculating flows 

By Y. Zang 1 , R. L. Street 1 AND J. R. Koseff 1 

The dynamic subgrid-scale model of Germano et al. (1991) is implemented in a 
finite volume formulation and applied to the simulation of turbulent flow in a three- 
dimensional lid-driven cavity at Reynolds number of 7500. The filtering operation 
is carried out in physical space, and the model coefficient is calculated locally. The 
computed mean and r.m.s. velocities as well as the Reynolds stress are compared 
with experimental data. It is shown that backscatter from small to large scales 
is necessary to sustain turbulent fluctuations. The model is being applied to the 
simulation of turbulent flows in stratified and rotating environment in complex 
geometries. 

1. Motivations and objectives 

Simulation of most engineering, environmental, or geophysical flows requires the 
use of turbulence model to represent unresolved small scale motions. Large eddy 
simulation (LES), which computes the motion of large spatial scales explicitly and 
models that of subgrid scales, proves to be an effective tool. During the last three 
decades, extensive effort has been put into developing subgrid-scale (SGS) turbu- 
lence models. Conventional SGS models such as the popular Smagorinsky (1963) 
model have several significant drawbacks. For example, they usually require model 
constants as input. These constants are not universal and have to be optimized to 
fit each type of flow. Although the Smagorinsky model has been successfully used in 
both homogeneous and wall-bounded turbulent flows (Bardina et al. 1983; Piomelli 
et al. 1988), it was found to be too dissipative in wall-bounded transitional flows 
(Piomelli et al. 1990), partly because the model does not account for the effect of 
energy backscatter from small to large scales. In addition, the Smagorinsky model 
gives a non-zero SGS stress in laminar flows. 

Recently, Germano et al. (1991) developed a new SGS model which calculates the 
model coefficient dynamically using the smallest resolved scales. This dynamic SGS 
model (hereafter referred to as DSM) has many desirable features such as requiring 
only one input parameter, exhibiting the correct asymptotic behavior near solid 
walls and in laminar flow, and being capable of accounting for energy backscatter. 
Most of the studies so far using the dynamic model has coupled it with spectral 
methods and computed flows with one or more homogeneous directions (Germano 
et. al. 1991; Moin et. al. 1991). 

The long-term goal of this project is to develop a numerical model to investigate 
environmental or geophysical flows in complex geometries. These flows are highly 
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inhomogeneous and may contain local regions of laminar, transitional, or turbulent 
flows of which successful simulations require a SGS model that can adjust itself to 
local flow dynamics - a quality which DSM possesses. The immediate objectives 
of the present work are to incorporate DSM into a finite volume Navier-Stokes 
code (Zang et al 1992) which is written in the general curvilinear coordinate 
system and to examine the model’s applicability to turbulent recirculating flows. 
We performed large eddy simulations of flow in a lid-driven cavity at a Reynolds 
number of 7500. The computed mean and fluctuating profiles were compared with 
experimental data. Local averaging of the model coefficient was investigated. The 
model is being applied to the simulation of stratified rotating upwelling flows on a 
slopping bottom. 


2. Accomplishments 

2.1. The finite -volume implementation of DSM 
The SGS Reynolds stress after applying the grid filtering operation represented 
by an overbar to the momentum equation is 

Tij = ITfUj - U t Uj ( 1 ) 

The Smagorinsky model is employed as the base model for as 

Tij - 6 fr kk = - 2 vrSij ( 2 ) 


where the eddy viscosity 

u t = CA 2 \S\ 


(3) 


A is the grid filter width, |S| = (2 S,jSij) 1/2 is the magnitude of the resolved large- 


scale strain rate tensor 


Sti 


iA 

2'dxj dxi 


(4) 


and C is the model coefficient which is to be computed dynamically using the 
procedure described below. We note that C is the square of the commonly used 
Smagorinsky constant Cs- 

Following Germano et al. (1991), we introduce a test-scale filter represented by a 

tilde whose length scale A is larger than A. The purpose of doing this is to utilize 
the information between the grid scale' and the test-scale filters to determine the 
characteristics of the SGS motion. By applying Germano’s algebraic identity (Ger- 
mano 1991) relating the subgrid-scale stress and the test-scale stress and employing 
a least-square closure (Lilly 1992), we obtain an explicit expression for the model 
coefficient C based on the local flowfield as 
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where 


Cij = UiUj — UiUj, (6a) 

Mij = a 2 |f|§ 0 - (66) 

a = A/A. (6c) 

In the above model, a is the only input parameter. 

The implementation of the SGS stress term l —drij/dxfl in the momentum equa- 
tions needs some consideration. First the trace of t X) is combined with the pressure 
gradient. Then, if the time marching scheme is explicit, Tjj can be calculated di- 
rectly from equations (2), (3) and (4). However, if an implicit or semi-implicit time 
marching scheme is used, it is advantageous to advance the SGS stress term (or at 
least a part of it) implicitly. This requires us to express SGS stress term in terms 
of of the velocity as 


drjj d_ du x . dvr du 3 

dxj dxj UT dxj dxj dxi 

In deriving the above equation, we have utilized the continuity equation. The first 
term on the right-hand-side of equation (7) can be combined with the viscous terms 
and be advanced in time implicitly. The second term, on the other hand, poses 
difficulty if treated implicitly because by doing so, the three momentum equations 
are coupled and the resulting matrix is more complicated. In the present work, this 
term is advanced explicitly as an extra source term. Numerical experiments showed 
that this treatment did not affect stability in the present case. 

The governing equations are transformed into the computational space and dis- 
cretized on a non-staggered grid using a finite volume approach (Zang et a 1. 1992). 
Cartesian velocity and volume fluxes are used as dependent variables. The solution 
method integrates the governing equations in time semi-implicitly by solving the 
viscous terms with the approximate factorization technique. The pressure Poisson 
equation is solved using a multigrid method. The overall accuracy is second-order 
in both space and time. The use of the non-staggered grid layout and Cartesian 
velocities as dependent variables makes the filtering operations very easy to carry 
out when the dynamic eddy viscosity is computed. 

Most studies to date using DSM have carried out filter operations in wave space. 
Recently, effort has been made to evaluate the behavior of different kinds of filters 
in physical space (Lund 1991). In the present formulation, we employ a top-hat 
filter in physical space with the trapezoidal rule for the test-scale filter. Because 
of the non-staggered-grid layout, all the variables that are used in the calculation 
of C are defined at the center of a control volume, and thus the filtering operation 
is the same for every variable. The length scale of the test filter is twice of that 
of the grid cell resulting in an a of 2. This value of a has been shown to be the 
optimal choice in the simulation of a turbulent channel flow (Germano et aL 1991). 


88 


Y. Zang, R . L. Street & J. R. Koseff 



FIGURE 1. Geometry and boundary condition of the lid-driven cavity flow. 

After the model coefficient C is computed from equation (4), it is averaged locally 
in space using the same scheme as that of the test filtering. 

When local averaging of C alone was employed, numerical instability was expe- 
rienced. This occurred because a large and negative C could lead to a negative 
eddy viscosity vt with a magnitude larger than the molecular viscosity, resulting 
in a negative total viscosity. Similar instability has been observed in previous sim- 
ulations using the dynamic model and a local model coefficient C(x,y, z,t) (Cabot 
1991). When there are one or more homogeneous directions, a plane or line aver- 
age in those directions could be used to stabilize the calculation (e.g., Germano et 
ai 1991). In the present wall-bounded flows, however, there is no homogeneous 
direction. Larger averaging volume and/or averaging in time could be employed to 
reduce the magnitude of negative C s (Piomelli 1991), but they still do not generally 
guarantee stability. To overcome this difficulty, we set the total viscosity vt + v 
equal to zero wherever it becomes negative. This cutoff eliminates the negative 
total viscosity that causes the numerical instability and at the same time allows 
energy to backscatter from small to large scales. The local averaging before the 
cutoff serves to spread the effect of large negative values of C if they exist. For the 
case computed in this work, the cutoff does not affect the flow significantly as the 
magnitude of the negative vt is larger than v only in a very localized region near 
the upper right corner of the cavity. 

2.2. LES of flow in a lid- driven cavity 

Lid-driven laminar flows in a cavity have long been used as a benchmark to 
test the accuracy of numerical methods (e.g., Ghia et al. 1982; Perng & Street 
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1989). Although extensive experimental measurements have been performed for 
transitional and turbulent flows in cavities (Koseff & Street 1984; Prasad & Koseff 
1989), to the authors’ knowledge, there have been no detailed direct or large-eddy 
simulations. Part of the difficulty lies in that in its temporal evolution the flow goes 
through the whole range of laminar, transitional, and turbulent stages. Moreover, 
even at the fully developed state, the flow is highly inhomogeneous and may have 
both laminar and turbulent regions. Successful simulation of this flow requires a 
SGS model which is capable of adjusting itself to the local flow dynamics. 

We computed the lid-driven flow in a three-dimensional cavity with aspect ratios 
B:D:H of 1:1:0. 5. A schematic of the cavity together with the notations are given in 
Figure 1. The Reynolds number based on the lid velocity Ub and the cavity length 
B is 7500. Past experiments (Koseff & Street 1984) have shown that at Reynolds 
numbers higher than about 6000, instabilities occur near the downstream eddy. As 
the Reynolds number increases, the flow becomes increasingly turbulent near walls, 
and at Reynolds numbers higher than 10,000, the flow near downstream eddy be- 
comes fully turbulent. Thus, the present flow is in the transitional regime. In the 
simulation, we use a 64x64x32 grid, which is non-uniform in the streamwise and 
vertical directions but is spanwise uniform. Initially the fluid is at rest. Computa- 
tion was first carried out in a half of the cavity with a 64x64x18 grid. A symmetry 
boundary condition was specified at the mid-plane. After the flow was fully de- 
veloped, the half flowfield was mirrored onto the whole cavity. Simulation then 
continued and statistics were collected after a relaxation time of three turnaround 
time-scale T a of the cavity. The value of T a , which was about 1 minute in this case, 
was estimated as the time for a particle at the edge of the top boundary layer to 
travel back to its starting position in the cavity. 

Figure 2 shows the computed mean center-line velocity profiles on the mid-plane. 
Measurements by Prasad and Koseff (1989) are also shown. The computed statistics 
were collected during a period of 5 T a , while the experimental data was averaged 
over periods of 5T a to 10T a . The computation slightly over-predicts the thickness of 
the boundary layer on the upstream wall but captures the thickness of the bottom 
boundary layer and the magnitude of the two peak velocities fairly accurately. The 
averaged relative error between the two data sets scaled by the lid velocity is less 
than 2 percent. 

The r.m.s. of fluctuating velocities \J < u ,2 >/Uq and \/ < v' 2 >/Ub and the 
Reynolds stress < tiV > /U B on the center-lines in the mid-plane are shown in 
Figure 3, where < • > denotes time averaging. We note that in the figure, these two 
quantities are increased by a scale of 10 and 500, respectively. The computation 
slightly over-predicts the r.m.s. of u r near the bottom boundary layer and under- 
predicts the r.m.s. of v* near the upstream wall. There is a minimum in the 
r.m.s. profile of u velocity in the bottom boundary layer which is not evident in the 
experimental data. The magnitude of the Reynolds stresses near the bottom and the 
upstream boundary layers is well represented. A maximum of the Reynolds stress is 
observed above the bottom boundary layer which is not shown in the experimental 
profile. 
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FIGURE 2. Mean centerline velocities on the mid-plane. Symbols are from Prasad 
and Koseff (1989). Lines are from the present computation. 

The computed profiles in Figures 2 and 3 are the large-scale quantities resolved 
by the grid, while the experimental data contains contribution from both the large 
and small scales. The contribution of the SGS motion to the Reynolds stress < 
uV > can be estimated using the time average of the SGS stress, < r 12 >. In 
the present case, the value of < t\ 2 > is at least an order of magnitude smaller 
than < v!v* > which is the large-scale contribution to the Reynolds stress. This 
indicates that the time-averaged statistics are well represented by the large-scale 
quantities. Two factors besides modeling error could have contributed to the above 
discrepancies. One is the experimental uncertainty, and the other is the effect of 
numerical resolution. However, at the present time, collecting statistics on a grid 
substantially finer than the one presently used is prohibitively expensive. 

Figure 4a shows the dynamically computed C on the mid-plane of the cavity, 
and 4b displays C on a plane close to one of the side-walls. On the mid-plane, 
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FIGURE 3. R.m.s. velocity and Reynolds stres s at t he centerlines on the mid- 
plane. Urms = 10\/ < u' 2 >/U B , V rm , = 10y/< V 2 >/U B , uv = 500 < u'v' > 
/Uq. Symbols are from Prasad and Koseff (1989). Lines are from the present 
computation. 

the range of C is from -0.12 to 0.1, while on the near- wall plane, the range is from 
-0.17 to 0.22. In the bulk region of the mid-plane, the magnitude of C is from 0.01 
to 0.02, which is comparable to the square of the commonly used value of 0.1 for 
the Smagorinsky constant. On the other hand, near the side wall (Figure 4b), C is 
small except near the corners and in the corner boundary layers. This is consistent 
with the expected behavior of C near a solid wall. It is interesting to notice that C 
is small near the moving top lid in both planes. 

There axe localized regions in Figure 4 where C is negative, which results in 
negative eddy viscosity representing energy backscatter from small to large scales. 
The ability of the present model to backscatter energy to large scales is important 
in sustaining turbulent fluctuations in the simulation. The time history of the 
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FIGURE 4a. Contours of computed C on the mid-plane. Dotted lines represent 
negative values. 



FIGURE 4b. Contours of computed C near one of the side walls. Dotted lines 
represent negative values. 


large-scale streamwise velocity u(0 near the peak of the bottom boundary layer on 
the center-line of the mid-plane is shown in Figure 5a, where the flow was in its 
fully developed state. Figure 5b gives tl(/) at the same location when no negative 
ur was allowed. The fluctuations in 5b were slowly damped out, indicating that 
backscatter from small to large scales is necessary to sustain turbulence. The low 
frequency oscillations with a period of about 1 minute in Figure 5a correspond 
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FIGURE 5a. Time history of u near the peak of the lower boundary layer on 
the vertical center-line in the mid-plane. Model with backscatter. Fully developed 
state. 



FIGURE 5b. Same as Figure 5a, Model without backscatter. 

to the spanwise meandering of TGL vortices. Examination of the flow structures 
in the spanwise direction shows the appearance of one pair of TGL vortices near 
the downstream wall which is meandering spanwisely. This is consistent with past 
experimental results (Prasad 1989). 

Previous simulations of cavity flows at lower Reynolds numbers have shown that 
the flows are essentially symmetric over the mid-plane (Perng &; Street 1989). This 
is true when the flow is laminar since there is no non-symmetric forcing and any 
non-symmetric small disturbances are damped. However, in turbulent flows, small 
disturbances could be amplified and result in asymmetry. In the present case, it 
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was found that when the half-cavity domain was used and the symmetry boundary 
condition was imposed at the mid-plane, both the streamwise fluctuating velocity 
< u' 2 > and the streamwise- vertical Reynolds stress < u'v ' > were unphysically 
large. On the other hand, the mean velocity profiles were barely changed. This 
was because the symmetry boundary condition eliminated the spanwise fluctuating 
velocity w at the mid-plane and restricted the meandering of the TGL vortices 
which were responsible for the momentum and energy transfer in the spanwise 
direction near the mid-plane. These high frequency fluctuations are averaged out 
in the mean profiles but make significant contribution to the r.m.s. quantities and 
Reynolds stress. 

3. Summaries and future plans 

A dynamic SGS model is coupled with a finite volume solution method and 
employed in the large eddy simulation of turbulent flow in a lid-driven cavity. Local 
averaging together with a cutoff is employed to obtain the model coefficient. The 
mean and fluctuating quantities were compared with experimental data and good 
agreement was achieved. It was shown that backseat ter is necessary to sustain 
fluctuations in the flow. 

The model is being applied to the simulation of upwelling flows of a stratified ro- 
tating fluid on a slopping bottom. Baroclinic instability was observed at the surface 
density front and intensive mixing occurred on both sides of the front. Preliminary 
results have shown that the computed wave speed and wave size compare favorably 
with the experimental and theoretical values. The characteristics of the instabilities 
and the subsequent breakdown of the front are being investigated. The model is 
also being employed to investigate flows in more complex geometries. 
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